Motion of Wavefronts in Semiconductor Superlattices 
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An analysis of wavefront motion in weakly coupled doped semiconductor superlattices is pre- 
sented. If a dimensionless doping is sufficiently large, the superlattice behaves as a discrete system 
presenting front propagation failure and the wavefronts can be described near the threshold currents 
Ji (i = 1, 2) at which they depin and move. The wavefront velocity scales with current as J — Ji | 3 . 
If the dimensionless doping is low enough, the superlattice behaves as a continuum system and 
wavefronts are essentially shock waves whose velocity obeys an equal area rule. 
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I. INTRODUCTION 



Semiconductor superlattices (SL) are unique nonlin- 
ear systems. Experimental evidence shows features of 
spatially discrete systems (multistability of roughly as 
many stationary states as SL periods due to formation of 
electric field domains Q) in certain regimes of nonlinear 
charge transport in SL. In other regimes (self-sustained 
periodic and chaotic oscillations of the current ||), SL 
behavior is more typical of continuous systems. These 
different properties may be described by discrete || or 
continuous balance equation models Q]. The continuum 
limit of discrete SL models has been used to understand 
important aspects of self-oscillations On the other 
hand, discrete models of SL share common characteris- 
tics with spatially discrete systems in other fields: most 
importantly, front propagation failure and front depin- 
ning @. 

Discrete systems describe physical reality in many dif- 
ferent fields: propagation of nerve impulses along mieli- 
nated fibers Jt],[| , pulse propagation through cardiac cells 
H , calcium release waves in living cells |J, sliding of 
charge density waves jl0| , superconductor Josephson ar- 
ray junctions pd| , arrays of coupled diode resonators Jl2| , 
motion of dislocations in crystals 1 13 and atoms adsorbed 
on a periodic substrate G3] . A distinctive feature of dis- 
crete systems (not shared by continuous ones) is the phe- 
nomenon of wavefront pinning: for values of a control 
parameter in a certain interval, wavefronts joining two 
different constant states fail to propagate Q . When the 
control parameter surpasses a threshold, the wavefront 
dcpins and starts moving |J^,|I^,|l|. The existence of 
such thresholds is thought to be an intrinsecally discrete 
fact, which is lost in continuum aproximations. The char- 
acterization of propagation failure and front depinning 
in discrete systems is thus an important problem, which 
is still poorly understood despite the numerous inroads 
made in the literature @,||l0]ji]| . 

Recently, several of us have proposed a theory of front 



depinning and propagation failure in discrete reaction- 
diffusion (RD) systems This theory will be extended 
here to describe motion of wavefronts near depinning 
thresholds in discrete drift-diffusion (DDD) models of SL 
p6[ . There are important differences between RD and SL 
systems. Let us consider a RD system depending on a 
parameter A, measuring how close we are to the contin- 
uum limit {A = 0), and an external parameter ('the field' 
F) such that wavefronts move for all values of F ^ in 
the continuum limit. In a generic RD system, there is 
always a pinning interval (wavefronts are stationary for 
F in the pinning interval) about F — no matter how 
close to the continuum limit we are. As we approach 
the continuum limit, the pinning interval shrinks to zero 
exponentially fast as A approaches zero: we then need 
exponential asymptotics (asymptotics beyond all orders) 
to describe what happens. In a SL, there are important 
differences: (i) the pinning interval disappears before we 
reach the continuum limit, and (ii) the continuum limit 
is described by hyperbolic equations and shock waves (no 
exponential asymptotics is needed). 

Let us consider an infinitely long, sufficiently doped SL 
under constant current conditions. The current density 
plays the role of external field in the SL models and a 
dimensionless doping v plays the same role as the pa- 
rameter A: the continuum limit is reached as v — > 0. For 
large values of v, we are in the strongly discrete limit and 
wavefronts are stationary (pinned) for currents on a cer- 
tain pinning interval. If the current is smaller than the 
lowest limit of the interval, J\{y), the wavefront moves 
downstream, following the electron flow. For currents 
larger than the upper limit of the interval, Ji [y] , wave- 
fronts move upstream, against the electron flow ||. As 
the doping decreases, first Ji and then J\ (for even lower 
doping) disappear. This means that fronts are either sta- 
tionary or move downstream below a first critical doping 
and they always move downstream below a second critical 
doping. Since the continuum limit corresponds to van- 
ishing dimensionless doping, stationary fronts and fronts 
moving upstream are features of discrete SL, which are 
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lost in the continuum limit. 

In contrast with the precipituous jumps and falls of 
charge and electric field in individual wells during wave- 
front motion, we will show that, for currents outside 
(Ji, J2), wavefronts are described by continuous profiles 
moving at a constant velocity. These profiles become 
sharper as the currents approach the threshold values J^; 
see Fig. |l|. Exactly at these values, wavefront profiles 
cease to be continuous: a number of jump discontinuities 
open up, which results in propagation failure. We de- 
scribe propagation failure by analyzing the behavior of 
several active wells which govern front motion. Our the- 
ory yields the front profile near threshold currents and a 
universal scaling of the front velocity for sufficiently large 
doping. Its performance worsens as the doping decreases 
and we approach the continuum limit. We will comple- 
ment our approach by obtaining the velocity of wave- 
fronts and their shape in this limit. Previously, these 
front profiles and velocities in the continuum limit were 
known only for models having no discrete diffusion [ jl7| , 
which is correct only for high enough voltage bias. For 
bias on the first plateau of the current-voltage charac- 
teristics, diffusivity cannot be ignored. Given that char- 
acterizing wavefronts at constant current is the key to 
analyze self-oscillations at constant voltage, the results 
in this paper could be useful to describe self-oscillations 
at low biases. This could be done in two limits: near 
the continuum limit !/<l, the results we have obtained 
for shock waves and monopoles can be used to charac- 
terize self-oscillations, either mediated by monopoles || 
or by dipoles |Q, For larger v, our present theory 
paves the way to understanding self-oscillations in DDD 
or more general SL models. 

The rest of the paper is as follows. We describe front 
propagation failure and front depinning for DDD SL 
models in Section ||. The continuum limit of the model 
is analyzed in Section HI. The resulting equations are 
those of the Kroemer model of the Gunn effect with zero 
diffusivity Wavefronts are shock wave solutions of this 
model with a rigid tail region either to the right or left of 
the shock. The shock velocity must be obtained from the 
discrete model, which results in an equal-area rule rather 
different from that of the Gunn effect. Corrections to this 
rule are important for describing the self-oscillations and 
will be given here. Section IV contains our conclusions. 



Lastly, the Appendices are devoted to different technical 
matters. 



II. PROPAGATION FAILURE AND FRONT 
DEPINNING 

In weakly coupled SL, (inter and intra-subband) scat- 
tering times are much shorter than well escape times, 
which in turn are much shorter than macroscopic times 
(period of self-oscillations). Then the dominant mecha- 
nism of vertical transport is sequential tunneling, only 



the first subband of each well is appreciably occupied 
and the tunneling current across barriers is stationary 
Q. Nonlinear phenomena seen in experiments can be 
described by discrete balance equations. 

We count the barrier separating the injecting contact 
from the first well of the SL as the zeroth barrier. Barri- 
ers and wells have widths d and w, respectively, so that 
I = d + w is the SL period. The ith SL period starts just 
before the ith barrier and ends just before the (i + l)st 
barrier. With this convention, the dependent variables of 
our model are minus the electric field averaged over the 
ith period, Fi, and the two-dimensional electron density 
at the ith well (concentrated in a plane normal to the 
growth direction, located at the end of the ith. well), rii. 
These variables obey the Poisson and charge continuity 
equations: 



Fi - 



drii 
~dt 



Ji-l- 



-N'd). 

Ji— >i+l • 



(1) 

(2) 



Here Np, e, e and eJ^i + i are the 2D doping density 
at the ith well, the average permittivity of the SL, mi- 
nus the electron charge and the tunneling current density 
across the ith barrier, respectively. We can differentiate 
(0) with respect to time and eliminate by using (0). 
The result can be written as a form of Ampere's law for 
the balance of current 



£ dF\ 
e dt 



Ji- 



= J(t). 



(3) 



Here eJ(t) is the total current density through the SL, 
equal for all SL periods, and e dFi/dt is the displacement 
current at the ith SL period. 

The tunneling current density eJi^+i is related to 
electric fields and electron densities by a constitutive re- 
lation, which should be derived from first principles; see 
Wacker in Ref. ||] and Rcf. |jq] . At sufficiently low 
or high temperatures, Ji—,i+\ has the following drift- 
diffusion form [fl6|, 



Ji — >i-\ 



niv(Fi) _ n l+1 - ni 



I 



P 



(4) 



where the drift velocity is an odd function of the field, 
v{— F) — —v(F), and the diffusion coefficient satisfies 
the relation 



D(-F) = v{F)l + D(F). 



(5) 



Typical forms of these coefficients (in nondimensional 
form; see below) are shown in Figure |^. To compare our 
theoretical results with numerical solutions of the model, 
it is better to use analytical approximations of these func- 
tions. Ours are given in Appendix |A|. The Ampere's law 
corresponding to the current (W) is 



e dFi niv(Fi) 



dt 



I 



-D(Fi 



n i+ i - Uj 
P 



= J(t) 



(6) 
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Equations (0) and @ form the DDD model stud- 
ied here. Notice that for high fields (in practice for all 
plateaus except the first one), D = 0, and we have a 
simpler discrete drift model with J^i+i = riiv(Fi)/l ||. 

To analyze the discrete drift-diffusion model, it is con- 
venient to render all equations dimcnsionless. Let v(F) 
reach its first positive maximum at (Fm,vm)- We adopt 
F M , N%, v M , v M l, eN%v M /l and eF M l/(eN%v M ) as 
the units of Ej, m, v(F), D(F), eJ and t, respec- 
tively. For the first plateau of the 9/4 SL of Ref. M, 
we find F M = 6.92 kV/cm, N% = 1.5 x 10 11 cm", 
vm = 156 cm/s, vm (d+ w) = 2.03 x 10~ 4 cm 2 /s and 
eN%v M /(d+w) = 2.88 A/cm 2 . The units of current and 
time are 0.326 mA and 2.76 ns, respectively. Then (|l|) 
and (^) become 

dE 

+ v{Ei) n, - D{Ei) (n;+i - m) = J, (7) 



Ei - Ej_! = i/ (nj - 1), 



(8) 



or, equivalcntly, 



dE 1 ,; ,„ , Ei — Ei_i ^ s Ej + i + Ej_i — 2E, t 



D(E Z ) 



J-v(Ei). (9) 



Here we have used the same symbols for dimensional and 
dimensionless quantities except for the electric field (F 
dimensional, E dimensionless). v — eN^/(e Fm) is the 
dimensionless doping parameter, which is about 3 for the 
first plateau of the 9/4 SL. 



A. Phase diagram for wavefronts on an infinite 
superlattice 

Here we shall consider an infinite SL under constant 
current bias J. Clearly, there are two stable spatially 
homogeneous stationary solutions, namely E^fJ) and 
E<- 3 \J), where v(E^) = J, E«(J) < E^{J) < 
E@\J). We are interested in nonuniform front states of 
the DDD model which satisfy E, -> E {1 \j) as i — y — oo 
and Ei -> E( 3 )(J) as i — * oo. These states are either 
stationary or time-dependent. In the second case, they 
are wavefronts moving with constant velocity c = c( J, u), 
such that Ei(t) = E(i — ct), i = 0,±1,... E(r) is a 
smooth profile which solves the following nonlinear eigen- 
value problem for c (measured in wells traversed per unit 
time) and E(r): 

dE , . E - E(t - 1) 
c — = W (E) - J + w(E) i '- 



dr 



-D(E) 



E(r+1)+E(r-1)-2E 



E(-oo) =EW(J), E(oo)^E^(J). 



(10) 
(11) 



By using a comparison principle [£(J, we can prove the 
existence of stationary fronts rigorously H . Outside the 



interval of current values in which there are stationary 
fronts, we can only prove that there are fronts moving 
to the right or the left |6). By using the comparison 
principle, it is possible to show that moving and sta- 
tionary fronts cannot exist simultaneously at the same 
value of the current . This result does not hold if the 
spatially discrete differential equation is second order in 
time: moving and stationary fronts have been reported 
to coexist in one such case Q . 

Proofs that moving fronts are traveling wavefronts 
have been given for different spatially discrete models 
without convection . In the case of the DDD SL 

model, we rely on numerical evidence: Notice that all 
noticeable steps in Fig. |l|, particularly (a) and (c), are of 
the same length. Fig. [| shows the wavefront at three dif- 
ferent times, and it clearly demonstrates that the front 
is a traveling wave moving with constant velocity as a 
whole. Our asymptotic construction of the wavefronts 
near the critical currents J\ and J2, which we explain 
below, exploits their traveling wave nature. 

Solving numerically (||), it can be shown that, after a 
short transient, a variety of initial conditions such that 
Ei -> EW(J) as % -► -00 and E l -> E^{J) as i — y 00 
evolve towards either a stationary or moving monopole. 
Figure 3 of Ref. [§ is a phase diagram showing the re- 
gions in the plane (J, v) where different fronts are stable. 
There are two important values of v, v\ < V2, such that: 

• For < v < v 1 and each fixed J £ (v m ,l), only 
traveling monopole fronts moving downstream (to 
the right) were observed. For v > u\, stationary 
monopoles were found. 

• For v\ < v < i>2, traveling fronts moving down- 
stream exist only if J g (v m , Ji(y)), where J\{v) £ 
(v m , 1) is a critical value of the current. J\{v) is a 
monotone decreasing function such that J\ — > v m 
as v — > 00; see Fig. [| If J e (Ji(v), 1), the stable 
solutions are steady fronts (stationary monopoles). 

• New solutions are observed for v > v^. As be- 
fore, there are traveling fronts moving downstream 
if J £ (v rn , Ji(is)), and stationary monopoles if 
J G (J\(v), J2{v)) 1 Ji(y) < 1 is a new critical 
current; see Fig. m The function Jzty) starts at 
^2(^2) = 1, decreases to a single minimum value, 
and then it increases towards 1 as v — > 00. For 
J%(y) < J <\, the stable solutions of ( p^ ) and ( |l3| ) 
are monopoles traveling upstream (to the left). As 
v increases, J\(y) and Ji(y) approach v m and 1, 
respectively. Thus stationary solutions are found 
for most values of J if v is large enough. 

Wavefront velocity as a function of current has been 
depicted in Figure g| for v = 3, which corresponds to the 
first plateau of the 9/4 SL. For larger v, the interval of J 
for which stationary solutions exist becomes wider again, 
trying to span the whole interval (v m , 1) as v — > 00. For 
very large i/, the velocities of downstream and upstream 
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moving monopoles become extremely small in absolute 
value. 



B. Pinning of wavefronts with a single active well 

At the critical currents, Ji{y) and J2(^), wavefronts 
moving downstream (to the right, following the electron 
flow, c > 0) for smaller J or upstream (to the left, against 
the electron flow, c < 0) for larger J fail to propagate. 
What happens is that the wavefront field profile E(t) 
becomes sharper as J approaches the critical currents. 
Exactly at Jk, k = 1,2, gaps open up in the wavefront 
profile which therefore loses continuity. The resulting 
field profile is a stationary front E^ = Ei(J, v): the wave- 
front is pinned for J\ < J < J2. The depinning transition 
(from stationary fronts to moving wavefronts) is techni- 
cally speaking a global saddle-node bifurcation. We shall 
study it first in the simplest case of large dimensionless 
doping ^, and then indicate what happens in the general 
case. 

For sufficiently large doping and J close to critical, the 
moving front is led by the behavior of a single well, which 
we will call the active well. If we examine the shape of a 
stationary front near the critical current, we observe that 
all wells are close either to E^ x \j) or E^\j) except for 
one well which drifts slowly and eventually jumps: the 
active well. Let us call E$ the electric field at the active 
well. Since all wells in the front perform the same mo- 
tion, we can reconstruct the profile E[i — ct) from the 
time evolution of Eo(t) = E(—ct). Before the active well 
jumps, Ei w £ (1) (J) for i < and E t w E (3 \j) for 
% > 0. Thus Eq. (||) becomes 

dE ° ~ 7 fj?> (W , E 0- E{1) 

a j^ v{ Eo)-v{Eo) 



-D(E ) 



+ £?(3) _ 2E 



(12) 



This equation has three stationary solutions for J\ < J < 
J2, two stable and one unstable, and only one stable sta- 
tionary solution otherwise. At the critical currents, two 
of these solutions coalesce forming a saddle-node. At low 
values of the current, the two coalescing solutions form 
a double zero corresponding to a maximum of the right 
side of ([l2]). For high currents, the two coalescing solu- 
tions form a double zero corresponding to a minimum of 
the right side of (12). The critical currents are such that 
the expansion of the right hand side of ([l2]) about the 
two coalescing stationary solutions, 



J - v(Eo) - v(E ) 



Eq - £?d) 



-D(E 



EW + £( 3 ) - 2E 



0. 



D' (EW + - 2E ) - 2D 
-v o (E -E^)-v o -isv' o = 0. 



(14) 



Here D' means D'(E ), etc. Equations (|l|) and @) 
yield approximations to Eq and the critical current J c 
(which is either J\ or J2). The results are depicted in 
Figs. ||and [| and show excellent agreement with those of 
numerical solutions of the model for v > 2. Our approx- 
imation performs less well for smaller i>, which indicates 
that more active wells are needed to improve it. 

Let us now construct the profile of the traveling wave- 
fronts after depinning, for J sligthly below J\ or slightly 
above J2 ■ Up to terms of order | J — J c \ , Equation ( |l2| ) 
becomes 



— w a (J 

dt K 



Jc 



Pf 2 



(15) 



for E Q (t) = Eq(J c ) + ip(t), as J — > J c . E (J C ) is the 
stationary solution of (Oh at J — J c . The coefficients a 
and are given by 



vq + Dq 



V Vn 



2v0 = D%(EW 



E {3) - 2E ) - AD'q - 2v' 
+v' Q \E^ -Eq-2v). 



(3 is negative if J c — J\ and positive if J c = J?,. 
has the (outer) solution 



(16) 

(17) 
Eq. @ 



<P® - (-l) ^^ ^ tan (V«/3(J- J fc ) (i - t )) 

(18) 

(k = 1,2), for J such that sign(J — J^) =sign/3. jl^ ) 
is very small most of the time, but it blows up when 
the argument of the tangent function approaches ±7r/2. 
Thus the outer approximation holds over a time inter- 
val (t — to) ~ ir/y/a/3 (J — Jfe). The reciprocal of this 
time interval yields an approximation for the wavefront 
velocity, 



\c(J,u)\ 



MR) 



(19) 



has zero linear term, 



In Figs. ^ and q, we compare this approximation with 
the numerically computed velocity for v = 3 and v = 20, 
respectively. 

When the solution begins to blow up, the outer solu- 
tion (|l^) is no longer a good approximation, for Eo(t) 
departs from the stationary value Eq(J c ). We must 
go back to ( [l2] ) and obtain an inner approximation to 
this equation. As J is close to J c and E (t) — E (J C ) 
is of order 1, we solve numerically (jlj) at J = J c 
with the matching condition that Eo(t) — Eq(J c ) ~ 
(!3) (-l) k 2/[nyp/[a(j-J c )] - 2|/3| (t - t )}, as (t - i ) ^ 
—00. This inner solution describes the jump of Eq to 
values close to E^ if J c = J\, or to values close to E^ 
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if J c = J2. During this jump, the motion of E forces the 
other points to move. Thus, for J c — J\, Ei(t) can be 
calculated by using the inner solution in (0) for Eg , with 
J = J c and E% « E^ 3 \ Similarly, for J c = J 2 , E-i(t) can 
be calculated by using the inner solution in (||) for Eq, 
with J = J c and i?_2 ~ -E^ 1 -*. A composite expansion 
p4| constructed with these inner and outer solutions is 
compared to the numerical solution of the model in Fig. 
[?]. Notice that we have reconstructed the traveling wave 
profiles E(i — ct) from the identity E (t) = E{—ct). 



C. Pinning of wavefronts with several active wells 

The previous asymptotic description of the depinning 
transition deteriorates as v decreases. What happens 
is that we need more than one active well to approxi- 
mate wavefront motion. Depinning is then described by 
a reduced system of more than one degree of freedom 
corresponding to active wells. There is a saddle-node 
bifurcation in this reduced system whose normal form is 
again ([l5]) with different coefficients. The jump of the ac- 
tive wells after blow up is found by solving the reduced 
system with a matching condition. How do we deter- 
mine the optimal number of active wells? For a given u, 
the eigenvector corresponding to the zero eigenvalue has 
a certain number of components that are of order one, 
whereas all others are very small. The number of com- 
ponents of normal size determines the optimal number 
of active wells: only one well if v is larger than 10, four 
if v = 3, etc. The eigenvector of the reduced system of 
equations for the active wells is a good approximation to 
the large components of the eigenvector corresponding 
to the complete system. As we approach the continuum 
limit, more and more points enter the reduced system 
of equations and the approach of Section III becomes a 
viable alternative to these methods. 

As before, the wavefront profile will obey E(i — ct) — 
Ei(t), where i = —L,...,M are the indices of the 
A = L + M + 1 active wells. Then their reduced dy- 
namics obeys (^|) for i = —L, . . . , M together with the 
approximations E_l-i ~ E^ x > and Em+i ~ E^ . We 
have to find approximate stationary solutions of this sys- 
tem, and then linearize it with respect to the appropriate 
one. At the approximate critical current J c , one of the 
eigenvalues of the linearized system becomes zero. Let us 
call U} and U l the A-dimensional (A = L + M + 1) left 
and right eigenvectors corresponding to the zero eigen- 
value of the coefficient matrix for the linearized system, 



chosen so that 



U}Ui 



1. Very close to J c , the 



electric field profile will be Ei(t) ~ Ei(J c , i/) + Ui<p(t) plus 
terms which decrease exponentially fast as time elapses. 
if obeys Eq. (R5h with the following coefficients 



a 



M 

E 



u, 



t . "0 + Do 



-L-l 



DmUI 



J M+l 



(20) 







M 



2v ^ 



[D'!U*{E i+1 +E i _ 1 -2E i ) 



+2DiUi(U i+1 + E/U - 2Ui) - 2vlUi(Ui - U^) 

+vZ(E^ 1 -E i -2v)U?]ul (21) 



Here we have set U-l-i 
E M+1 = E&, v, = v(E, 
All quantities are evaluated at J 



= U M +i = 0, - 
,„ v'i = v'(Ei) Di = D(Ei), etc. 



J c . Clearly these for- 



mulas become (|T|) and @ if L = M = Q, U = f/ f = 1. 

Obviously, the solution and interpretation of ( |l5| ) is the 
same as before, and we only have to change the matching 
condition and the description of the jump of the active 
wells. Increasing the number of active wells improves 
appreciably the approximation of the wavefront velocity; 
see Fig. § The jump of the active wells is described 
by the solution of (||) for i = —L, . . . , M, with J = J c , 
Ei(t) = Ei(J c , v) + Ui^p(t) 1 and the same matching con- 
dition as before. Fig. || depicts a comparison of the com- 
posite approximation and the numerically calculated field 
profile for v = 3. Notice the improvement with respect 
to the single active well approximation. 

Our analysis shows that the transition from moving to 
stationary fronts involves loss of continuity of the wave- 
front profile. For the reduced system of active wells, the 
transition is a global saddle-node bifurcation: two sta- 
tionary solutions coalesce at the critical current and (for 
J < Jiiy) or J > J2{v)), a traveling wavefront appears. 
Near the critical currents, the field profile of this front 
is very sharp: it resembles the typical discrete steps of 
the stationary profile on long intervals of time (scaling as 
J— Jt I ~ 2 ) , followed by abrupt transitions between steps. 
Although we do not have a proof, it is natural to conjec- 
ture that the depinning transition for the complete SL 
model is of the same type as that for the reduced system 
of active wells. 

Remark. Notice that our theory of front depinning 
applies with minor modifications to discrete models 
with more complicated (than (Q)) constitutive relations, 
Ji^i+i = J{Fi,ni,n i+ i) (cf Wacker's paper in Ref. ||]) 
or Ji^i+i = J{F i ^ 1 ,F l ^\ +1 ,n i ,n i+1 ) (cf the complete 
discrete model in Ref. |l6| where potential drops at the 
barriers, Vi, are used instead of electric fields averaged 
over one SL period, Fi). The key point is that only a 
finite number of wells are active in wavefront depinning 
for sufficiently large v. 



III. CONTINUUM LIMIT 

The continuum limit of the DDD model is useful to 
understand self-sustained oscillations of the current and 
wavefront motion. It consists of v — > 0, i —* 00, with 
vi = x <E [0, Nu], Nu 3> 1. In this limit, (|!j) becomes 



(22) 
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up to terms of order v. Equation ( |22| ) corresponds to 
the hyperbolic limit of the well-known Kroemer model 
of the Gunn effect ]l8| . With constant J, shock waves 
are solutions of these equations [p5[ . These waves are 
related to wavefronts, and their speed can be calcu- 
lated explicitly. Let V(E + ,E^) denote the speed of 
a shock wave such that E becomes E_ (resp. E + ) to 
the left (resp. right) of the shock. Inside the shock 
wave, we should use the discrete model with n, ^> 1, 
dEi/dt V (Ei - E^-Cj/v > 1, and J = 0(1). No- 
tice that we have rescaled the wavefront velocity so that 
V = cv is the correct velocity for the (rescaled) contin- 
uum profile E(x - Vt) = E(y (i - ct)). Then 



i/n i+ i 



y-v^ Ei — Ei-i \ - D{E.. 



v(Ei) + D(Ei) ^ v(Ei) + D(Ei) 

~ Ej-i D(Ei) [E i+X - Ei 



v(Ei) + D(Ei) 
This expression yields 



v{Ei) + D(Ei) 



V 



v v{EA (E t+1 -E t ) 
^ v( 



v{E z ) + D{Ei) 

Ej-Ej-! 
iEi) + D(Ei) 



which becomes 



f E+ V (E) , F 
JE- v(E)+D(E) w±j 



V(E +1 E_)- 

V + ' ; fE+ _ dE 

JE_ v 



(23) 



{E) + D(E) 



in the continuum limit. This expression is equivalent to 
the following weighted equal-area-rule 



? + v(E) - V(E+,E. 
E _ v(E)+D(E) 



dE = 0. 



(24) 



For D = 0, this formula reduces to that derived for the 
discrete drift model in Ref. |17|. An expression includ- 
ing the next correction to this formula is given in Ap- 
pendix |b[ There is only one value of J, J* , such that 
V = J with E- = £W(J) and E+ = E^(J). For 
J E (v m , J*), a wavefront joining E^ (J) to E^ (J) con- 
sists of a shock wave having E + = E^(J), and E_ such 
that V(EW(J),E-) = v(E_). Furthermore, to the left 
of the shock wave, there is a tail region moving rigidly 
with the shock and such that 



dE 

[ v (E)~V] — = J-v(E), 



(25) 



for negative £ = x — Vt, and E{— oo) = E^'(J), 
E(0) = E_. This whole structure (shock and tail re- 
gion) is called a monopole with left tail ||. Similarly, for 
J £ (J*, 1), a wavefront joining E^{J) to E^(J) be- 
comes a monopole with right tail. This monopole consists 
of a shock wave having E- = E^\j), and E+ such that 
V(E+,E^(J)) = v(E + ), and a tail region satisfying M) 



for positive £, with the boundary conditions E(0) = E + 
and E(oo) = E^(J) g. In conclusion, the wavefront 
velocity as a function of J is determined by the following 
equations: 

c(j) = y(£ (3) (J),£-), with 

v(E-) = V(E {3 \ J), E-), if v m <J<J*, (26) 

C(J) = F(£ + ,£: (1) (J)), with 
v(E+) = V{E+,E {1) (J)), if J*<J<1. (27) 

Notice that this C( J) is the limiting value of the rescaled 
wavefront velocity, c( J, v)v as v —> 0+. Wc have 
compared the continuum approximation of the wave- 
front velocity (in wells traversed per unit time, i.e., 
c(J,f) = C(J)/v, not rescaled) with the numerical solu- 
tion of the model for v — 0.01 in Fig. ^. The equal-area 
rule result corresponds to (|25| ) , J26| ) and ( p7|) and its max- 
imum difference with the numerical solution is about 17.6 
%. Notice that the corrections in Appendix |^ improve 
sign ificantly the result: the corrected equal-area result 
(B2) differs at most 3% from the numerical solution. 



Remark. A different strategy to derive a formula for the 
shock velocity V(E + , EJ) could be to retain the first or- 
der corrections to fl22]) and then use well-known proce- 
dures for partial differential equations. The first-order 
correction to (E2I) is 



dE ..dE 



D(E) 



v(E) 



-v(E) 
d 2 E 



dx 2 



(28) 



in which terms of order v 2 are ignored. As v — > 0, this 
equation has shock waves whose velocity obeys the equal- 
area rule 



E+ 



v(E)-V(E + ,E_) 



D(E) 



v(E) 
2 



dE = 0. 



(29) 



instead of (|24|). Notice that this rule gives the same re- 
sult as ( p4| ) for D = 0. What is wrong in this argument? 
Notice that we have to use the rescaled moving variable 
£ = (x — Vt)/v to derive (||j|). After rescaling, both con- 
vective and diffusive terms in ( p8| ) are the largest ones, 
of order v^ 1 ^> 1. But so are all the terms in the Tay- 
lor series {E, - E^jv = dE/dx - {v/2)d 2 E/dx 2 + 
(f 2 /6) d 3 E/dx 3 + . . ., which was used to derive (|28| ) in 
the first place. Thus this derivation ignores infinitely 
many relevant terms in the Taylor expansions of rii and 
and it yields the incorrect formula (p9|) as a result. 



IV. CONCLUSIONS 

We have analyzed the motion of wavefronts in discrete 
drift-diffusion models of nonlinear charge transport in su- 
perlattices. Moving wavefronts are profiles of the electric 
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field traveling rigidly with constant velocity. Propagation 
failure of the fronts occurs because the field profile loses 
continuity at the critical currents and becomes pinned at 
discrete sites. We have characterized propagation failure 
of wavefronts and, conversely, front depinning by sin- 
gular perturbation methods. These methods are based 
upon the fact that only a few active wells characterize 
wavefront motion for dimensionless doping v sufficiently 
large. In the continuum limit, as v tends to zero, more 
and more wells become active and a different approx- 
imation makes sense. In this limit, discrete equations 
turn into differential equations, and wavefronts turn into 
monopoles which are shock waves with a rigid tail region 
to their left or right. We have derived the equal area rule 
for such shock waves, and also its leading correction. Our 
different asymptotic theories perform well in their respec- 
tive domains of validity and approximate the motion of a 
wavefront at constant current on an infinite superlattice. 
Understanding this motion is the key to understanding 
more sophisticated phenomena occuring in semiconduc- 
tor superlattices under voltage bias conditions 
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2 ^ 



[f(Ei) — f(Ei_i)](Ei — Ei-i) 



J(E l _ 1 )(E l -E^ 1 ) 



fiE^m-E^). (Bl) 



We can repeat our derivation in Section [II step by step, 
but keeping now the correction terms (Bl) when sub- 
stituting integrals instead of Riemann sums. The result 
is 



rE+ y(E) 

v(e + ,e_) = Je - e v{e)+d{e] 

li 



dE + S 



N 



dE 



s 



N 



Sn — 



2 ^ 

IE 



D, 



£_ v(E)+D(E) 

A-i 



■Di 



Di 



Vi-i + A-i 

Vj-l 

Vi-i + A-i 
1 

vt-i + A-i 



f S D 

(Ei - Ei-i) 



(B2) 



{Ei-Ek-x), (B3) 



{Ei-Ei-i). (B4) 



In the numerator of (B2), we have ignored the term 
i/y^fa — J)/(vi + Di), which vanishes in the continuum 
limit. Notice that vi-\ ~ Uj — v[ (Ei — a similar 

relation for A-i and that v\ < (typically) and D[ < 
suggest that Sn < and Sd > 0. Thus we expect that 
the corrected shock velocity be lower than the leading 
order (Ea). 



APPENDIX A: DRIFT VELOCITY AND 
DIFFUSIVITY COEFFICIENTS 



To compare our theoretical results with those of nu- 
merical solutions of the model, it is better to use an- 
alytical approximations of drift velocity and diffusivity 
instead of using numerical data for them. We use the 
following dimensionless coefficients which approximately 
fit those used in Ref. 101: 



v(E) 



0.2684 E 



0.25 + (0.9862£;-0.85) 2 
°J™ 3-0.004, (A!) 



0.16+ (0.9862£- 15) 
D(E) = 2e- E2 . 

The velocity has a local maximum at E = 1, v = 1. 



(A2) 



APPENDIX B: FIRST ORDER CORRECTION TO 
EQUAL- AREA RULE 

A more careful derivation of the equal-area rule uses 
the trapezoid rule for Riemann integrals: 



f(E)dE^Y.f( E ^( E *- E *-i) 
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FIG. 1. Sharpening of wavefront profiles as the dimension- 
less current J approaches its critical values for v = 3. (a) 
J ~ Ji, (b) J > J 2 , (c) J « Ji, (d) J < Ji. 




FIG. 2. Electron drift velocity v(E) and diffusion coeffi- 
cient D(E) as functions of the electric field in nondimensional 
units. 
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FIG. 3. Electric field profile of a wavefront at three con- 
secutive wells. We have depicted the field at each well as a 
function of time, which illustrates the motion of the front at 
constant velocity. 



FIG. 5. Wavefront velocity as a function of current den- 
sity for v — 3. We have compared the numerically measured 
velocity to the results of our theory with one or four active 
wells. 




V v FIG. 6. Wavefront velocity as a function of current density 

for v = 20. We have compared the numerically measured 
FIG. 4. Upper and lower critical currents as functions of velocity to the results of our theory with one active well, 
the dimensionless v. 
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FIG. 7. Wavefront profiles near J = J2 for v = 20. We 
compare the results of matched asymptotic expansions with 
one active well and the numerical solution of the model. 
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FIG. 8. Wavefront profiles near J = J 2 for v = 3. We 
compare the numerical solution of the model with the result 
of matched asymptotic expansions with 5 active wells. Notice 
that the largest error source is the theoretical estimate of the 
time period (an error of 200 in a 3,000 period). 
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FIG. 9. Comparison of the equal-area rule (leading order) 
and corrected equal-area rule (including first order correc- 
tions) approximations to the wavefront velocity with the nu- 
merical solution of the model for v = 0.01. 
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